Resonating Valence Bond Wave Functions for Strongly Frustrated Spin Systems 
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The Resonating Valence Bond (RVB) theory for two-dimensional quantum antiferromagnets is 
shown to be the correct paradigm for large enough "quantum frustration" . This scenario, proposed 
long time ago but never confirmed by microscopic calculations, is very strongly supported by a new 
type of variational wave function, which is extremely close to the exact ground state of the J1 — J2 
Heisenberg model for 0.4 < J2/J1 ^ 0.5. This wave function is proposed to represent the generic 
spin-half RVB ground state in spin liquids. 
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The question whether a frustrated spin-half system is 
well described by a spin-liquid ground state (GS) - with 
no type of crystalline order - 25 years after the first pro- 
posal [|j is still controversial, mainly because of the lack 
of reliable analytical or numerical solutions of model sys- 
tems. For unfrustrated or weakly frustrated quantum 
antiferromagnets a deep understanding of the nature of 
the GS together with a quantitative description of the or- 
dered phase is obtained by including Gaussian quantum 
fluctuations over a classical Neel state. pLpf For sizeable 
frustration, instead, this description is known to break 
down. However, the short-range RVB state M does not 
prove a good starting point for the description of frus- 
trated models; it rather turns out to be the exact GS of 
ad hoc Hamiltonians. 0-CT 

As a prototype of a realistic frustrated two-dimensional 
system, which has been recently realized experimentally 
in Li2VOSi04 compounds, W we investigate the spin- 
half Heisenberg model with nearest ( J\) and next-nearest 
neighbor ( J 2 ) superexchange couplings: 



W = Ji5^Si-S,- + J a Y. 
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on an N— site square lattice with periodic boundary con- 
ditions. In the (J 2 = 0) unfrustrated case, it is well 
established that the GS of the Heisenberg Hamiltonian 
has Neel long-range order, with a sizable value of the 
antiferromagnetic order parameter, flq] However, varia- 
tional studies H have shown that disordered, long-range 
RVB states have energies very close to the exact one. It 
is therefore natural to imagine that by turning on the 
next-nearest neighbor interaction J 2 , the combined ef- 
fect of frustration and zero-point motion may eventually 
melt antiferromagnetism and stabilize a non-magnetic 
GS of purely quantum-mechanical nature. Indeed, for 
0.4 < J2/J1 ^S 0.6 there is a general consensus on the 
disappearance of the Neel order towards a state whose 
nature is still much debated. ]10| 

In a seminal paper, pT| Anderson proposed that a 
physically transparent description of a RVB state can 



be obtained in fermionic representation by starting from 
a BCS-type pairing wave function (WF), of the form 



IV'bcs) = exp 
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This WF is the GS of the well-known BCS Hamilto- 
nian with a given (real) gap function A& = A_fc pro- 
vided the Fourier transform //. of the pairing function, 
fij, satisfies: f k = A k /(e k + ^Jt\ + A 2 .), where e k = 
— 2[cosk x + cos k y ] is the free-electron dispersion. The 
non-trivial character of this WF emerges when we re- 
strict to the subspace of fixed number of electrons (equal 
to the number of sites) and enforce Gutzwiller projec- 
tion onto the subspace with no double occupancies: sin- 
glet pairs do not overlap in real space and this WF can 
be described by a superposition of Valence Bond (VB) 
states. pT]-[ll| Though this WF has been already studied 
for the pure Heisenberg model by several authors []l2|jl3| 
for Afe oc (cos k x — cos k y ), here we show that this type of 
RVB state represents an extremely accurate variational 
ansatz for frustrated systems. 

A definite symmetry is guaranteed to the projected 
BCS (p-BCS) state provided the gap function A^ trans- 
forms according to a one dimensional representation of 
the spatial symmetry group. A careful analysis |1J], sim- 
ilar to the one carried out in fllal, shows that the odd 



component of the gap function A k 



^k+(w,Tr) may 



have spatial symmetries different from those of the even 



component Aj 



k fe+(7r,7r)- 



We anticipate that the best 



variational energy is obtained when the former has d x 2_ y 2 
symmetry, whereas the latter either vanishes or it has d xy 
symmetry. In frustrated models, it is important to con- 
sider this generalization of the originally proposed WF, 
|ll| , p2]| because only in this way it is possible to reproduce 
correctly the phases of the actual GS configurations. In 
the unfrustrated case it is well known that such phases 
are determined by the so-called Marshall-sign: on each 
real space configuration \x) , the sign of the GS wave func- 
tion is determined only by the number of spin down in 



one of the two sublattices. This feature, rigorously valid 
for J2 = 0, turns out to be a very robust property for 
weak frustration (J2/J1 ^0.3). lla] However, it is clearly 
violated when frustration plays an important role. 




0.3 
J 2 /J, 

FIG. 1. Average sign, accuracy of the GS energy, and 
overlap between the GS and the p-BCS state (full dots) as a 
function of J2/J1, for N = 6 x 6. Empty dots are the Marshall 
sign (top panel), and the energy accuracy of a Neel ordered 
spin- wave WF M (middle panel). Lines are guides for the eye 
and the shaded region indicates the location of the expected 
transition point to the non-magnetic phase. 

In order to determine the best variational WF of this 
form we have used a recently developed quantum Monte 
Carlo (QMC) technique jl6| that allows to optimize a 
large number of variational parameters with modest com- 
putational effort. We first consider the largest square 
cluster N — 6 x 6 where the exact GS can be numer- 
ically determined by exact diagonalization (ED). In or- 
der to show the quality of the present WF when frustra- 
tion (J2/J1) is increased, we have computed the vari- 
ational energy, the overlap of the p-BCS wave func- 
tion with the exact GS, \ipo), and the average sign, 
defined for a generic variational state \ipv) as (S) = 
EJ(2#v}| 2 Sgn[(2# v }<2#o}] . The Marshall sign (i.e., 
(S) = 1 for J 2 = 0) is obtained using the p-BCS wave 
function, with only the d x 2_ y 2 component. However, for 
J2/J1 ^ 0.4, the phases of the WF are considerably af- 
fected by the strong frustration and only when a sizable 
d xy component is stabilized at the variational level, this 
property can be correctly reproduced, as clearly shown 
in Fig. [y. Remarkably, as shown in the same figure, this 
kind of WF is not only in qualitative agreement with the 
exact solution but it is also impressively accurate in the 
region J2/J1 ~ 0.45 ±0.05 of large frustration, where the 
overlap of the variational WF is improved by more than 
an order of magnitude with respect to the J2 = case. 
This fact implies that the GS in the strongly frustrated 
regime is almost exactly reproduced by a RVB wave func- 
tion, at least on clusters of this size. 



Further indications of the changes in the nature of the 
GS occurring by increasing the frustration ratio can be 
found in the ordering of states with different quantum 
numbers in the energy spectrum. JlTj This information 
can be easily accessed by ED, which has been performed 
on the 6x6 cluster for three representative values of the 
frustration ratio: J2/J1 — 0.2; 0.5; 0.8. In the case of 
Neel order the two lowest states of the finite-size spec- 
trum are a total-symmetric singlet and a triplet of mo- 
mentum (tt, 7t). This phase is clearly expected to occur 
for sufficiently small J2/J1. Analogously, in the large 
J2/J1 limit, the two sublattices decouple and a collinear 
state characterized by ferromagnetic stripes, staggered 
along the direction orthogonal to the stripe, is believed 
to prevail. |L8| In this case, the symmetry breaking im- 
plies that four classes of states with different spatial sym- 
metries become degenerate: the lowest representatives 
of these families are an s-wave and a d-wave singlet at 
zero momentum and two triplets at momenta (0, tt) and 
(•7T,0). Therefore the transition between these two or- 
dered phases implies that, by increasing J2/J1, the (tt, n) 
triplet should acquire a gap while the d-wave singlet and 
the (0, 7r) and (tt, 0) triplets should collapse onto the GS. 
The low-energy spectrum is shown in Fig. 0, suggesting 
that the reshuffling of the lowest energy levels in the sys- 
tem occurs at least in two steps: first the triplet levels 
lift, leaving room for a non-magnetic phase with a finite 
triplet gap and then the d-wave singlet collapses. 
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FIG. 2. Lowest energy states with given quantum num- 
bers referenced to the GS energy: Singlet zero momentum 
s-wave (empty dot) and d-wave (cross), singlet at momentum 
(0, tt) (empty square), singlet at momentum (tt, tt) (empty tri- 
angle), triplets at momentum (n,it) s-wave (full dot) and at 
momentum (0,7r) (full square). 

In principle either a homogeneous spin liquid or a VB 
crystal with some broken spatial symmetry is compati- 
ble with a triplet gap in the excitation spectrum. Indeed, 
there are several proposals that the GS for J2/J1 ~ 0.5 
may be spontaneously dimerized |1C|]19|J20|| with broken 
translation-rotation symmetry (columnar VB state) as 
it is likely to happen in a generalization of the present 
model. M Alternatively, a plaquette VB state IP1 



with broken translational symmetry but preserving rota- 
tional symmetry may be stabilized. Note that a plaque- 
tte VB state would imply the degeneracy of four singlet 
states at momenta (0,0), {ir, 0), (0, ii) and (7r,7r), while 
a columnar VB state |l(||lj|^(| would result from the 
mixing of four singlets with different quantum numbers: 
two translationally invariant s-wave and <i-wave states 
and the two singlets at momenta (it, 0) and (0, it). The 
plaquette scenario recently proposed (2^^] turns out to 
be very unlikely in this model because the lowest sin- 
glet of momentum (7r, it) lies very high in energy at all 
the couplings we have investigated. The presence of a d- 
wave singlet in the singlet-triplet gap, instead, has been 
also evidenced in the regime of strong frustration of the 
J1 — J2 model on the 1/5-depleted square lattice, where 
the ground state is believed to be spin liquid. [£3| On the 
other hand on the basis of these ED data, it is clearly 
impossible to establish whether the GS is dimerized or 
disordered. 
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FIG. 3. Upper panels: energy vs. decreasing variance 

for the p-BCS wave function with zero, one and two Lanczos 

iterations. ED and QMC g (full triangles), are shown for 

comparison. Lower panel: variance vs. J2/J1. In all the 

plots: TV = 6 x 6 (dots) and, TV = 10 X 10 (squares). 

In order to clarify this issue we have extended the cal- 
culation to much larger system size, by also employing a 
few Lanczos iterations over the starting variational WF. 
The stochastically implemented Lanczos technique is a 
new QMC method with very good convergence proper- 
ties when the initial WF well represents the actual GS. 
Jig ] This accuracy can be confirmed a priori even on large 
size, by studying the variance a 2 = ((ft 2 ) — {'H} 2 )/J 2 N 
of the energy, the variance being smaller (zero) for a bet- 
ter (exact) calculation. As shown in Fig. |3[ also in the 



10 x 10 cluster the variance as a function of J2/J1 behaves 
similarly to the 6x6 case, strongly suggesting that the ex- 
ceptional accuracy of the p-BCS wave function does not 
decrease for larger sizes. This is also confirmed by the ex- 
tremely well-behaved approach to the zero- variance limit 
with few Lanczos iterations (shown in the same figure), 
leading to an almost exact estimate of the GS energy 
even for J 2 = 0, when the accuracy is the lowest. 
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FIG. 4. S(q) for (from the lower to the upper curve) 
TV = 6x6, 10 x 10, and 14 x 14. Variational esti- 
mate (empty triangles), with one Lanczos iteration (empty 
squares), with two Lanczos iterations (empty dots). Full dots: 
variance-extrapolated values of S(n,n); large empty circles: 
ED results. Inset: size-scaling of the variance-extrapolated or- 
der parameter squared for J2/J1 = 0.5 (full dots) and J2 = 
(empty dots, Ref. uM)- The full triangle is the thermody- 
namic value for J2 = taken from Ref. M. 

Of course, the accuracy in the energy does not neces- 
sarily guarantee a corresponding accuracy in correlation 
functions. However, in the strongly frustrated regime, 
our approach is particularly reliable since the gap to the 
first excitation belonging to the same subspace of our 
best WF (with two Lanczos iterations) is bounded in 
all the most plausible cases (columnar, plaquette, non- 
degenerate singlet RVB) by the triplet gap (> O.Ui, 
|I^,g2[) and therefore is much larger than the estimated 
error in the total energy (~ O.OUi, see top panel of 
Fig. g). Indeed, as shown in Fig. [|, the comparison of 
the magnetic structure factor S(q) — (S q ■ S_ 9 ), with the 
exact result gives a clear indication that correlation func- 
tions obtained by the variational approach are essentially 
exact, indicating also the absence of long-range Neel or- 
der. p2| This fact is particularly evident because for all 
the lattice sizes considered, S(tt,tt) is slightly depressed 
by few Lanczos iterations, meaning that the exact value 
of the magnetic structure factor is bounded by the one of 
the p-BCS wave function, with sizable antiferromagnetic 
correlations but with no antiferromagnetic long-range or- 
der. Remarkably, correlation functions are smoothly de- 
pending on the energy variance, so that an estimate of the 



magnetic order parameter within ~ 10% can be achieved 
also for Ji = 0, where our singlet WF is not particularly 
accurate and the spectrum is gapless (see the inset of 
Fig. |). 

In order to investigate the existence of long-range 
dimer-like correlations, as in the columnar or the plaque- 
tte VB state, we have calculated the climcr-dimer corre- 
lation functions, A*j = (SfSJSfSf) - (S*S*)(S^). In 
presence of a broken spatial symmetry, the latter should 
converge to a finite value for large distance. This is 
clearly ruled out by our results, shown in Fig. H, with 
a very robust indication of the liquid character of the 
GS for J2/J1 — 0.5, which is correctly described by our 
variational approach. 
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Dimer-dimer correlation functions A t '. obtained 



FIG. 5 
by keeping fixed the position of the left-most bond (i,j) (dou- 
ble stick) and moving the bond (k, I) (single stick) along the 
indicated patterns, d is the Manhattan distance. 6x6 (left), 
10 x 10 (right); symbols as in Fig. 0. 

A totally symmetric spin-liquid solution proposed for 
this model in Ref. pj was actually rather unexpected af- 
ter the work of Read and Sachdev, pfl] providing argu- 
ments in favor of spontaneous dimerization. This conclu- 
sion was supported by series expansion po|Jl9[ and QMC 
studies included the one done by two of us. |22| It is clear 
however that it is very hard to reproduce a fully sym- 
metric spin liquid GS, with any technique, numerical or 
analytical, based on reference states explicitly breaking 
some lattice symmetry. |gj] We do not know whether nu- 
merical methods and/or series expansions can ever solve 
this controversial issue. However, at the time being, we 
can safely state that in order to have something different 
from a spin- liquid GS, one has necessarily to improve the 
quality of our spin- liquid variational WF, e.g., with small 
symmetry breaking terms; a possibility that we have at- 
tempted (small dimerizations or plaquette-like perturba- 
tions) without success. Indeed, this seems a very difficult 



task due to the tiny energy range left (~ 10 -3 J\ per site) 
by our variational WF. 

In conclusion, the spin-liquid RVB ground state, orig- 
inally proposed to explain high- Temperature supercon- 
ductivity, is indeed very plausible for strongly frustrated 
spin systems. We expect that the p-BCS Resonating Va- 
lence Bond wave function represents the generic varia- 
tional state for a spin-half spin liquid, once the pairing 
function fij is exhaustively parameterized according to 
the symmetries of the Hamiltonian. In particular, the 
p-BCS wave function can be easily extended to the case 
of topologically frustrated lattices - like the kagome or 
the pyrochlore lattices - as well as to frustrated models 
on the 1/5-depleted square lattice. 
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